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ABSTRACT 


The  CORDIC  (Coordinate  Rotation  Digital  Computer)  algorithm"' 

computes  certain  functions  such  as  the  sine,  cosine,  and  using  only 

additions  and  bit  shifting  operations. 

We  have  implemented  an  integer  math  CORDIC  algorithm  on  a  high 
speed  RISC  processor.  During  the  course  of  this  work,  we  identified  a 

convergence  problem  with  the  Vx^T^  CORDIC.  A  solution  to  this  problem  is 
presented  along  with  an  overview  of  this  algorithm. 
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I.  INTRODUCTION 

The  CORDIC  algorithm^  utilizes  a  series  of  rotations  on  a  two 
dimensional  vector  to  compute  the  following:  sin(z),  cos(z),  arc  tan(y/x), 

and  In  j^s  generalized  version  it  has  also  been  shown  to  have  the 

capability  of  performing  multiplication  and  division,  as  well  as  computing 

hyperbolic  functions,  and  Vx2-y2 

CORDIC  has  found  its  way  into  desk  calculators,  specifically,  the  HP- 
9100  series2;  moreover,  it  has  proven  useful  in  calculating  the  Fourier 
Transforms,  and  also  the  singular  values  of  a  matrix^.  The  algorithm  can  be 
implemented  either  in  software  or  on  a  single  digital  IC^- 

We  first  discuss  the  CORDIC  algorithm,  and  then  present  a  problem  we 
encountered  in  its  use.  Since  our  project  involves  real  time  control  and 
requires  an  extremely  small  computer,  we  are  using  integer  math  in  an  RTX 
2000  processor^  programmed  in  its  native  FORTH  language.  A  problem 

arose  in  the  evaluation  of  using  CORDIC.  We  characterize  the 

problem  and  present  our  solution. 


II.  THEORY 

The  main  working  equations  of  the  CORDIC  algorithm  can  be  related  to 
the  orthogonal  transformation  equations  used  to  rotate  a  two  dimensional 
vector.  Let  us  assume  our  original  vector  R  has  components  x  and  y.  The 
transformation  equations  which  rotate  this  vector  through  a  positive 
clockwise  angle  5  are  ; 

(1 )  x'  =  X  cos  (6)  +  y  sin(5) 

(2)  y’  =  -X  sin(5)  +  y  cos(5) 
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Figure  1 :  Orthogonal  Rotation 


Since  the  polar  coordinate  0  of  a  vector  is  normally  defined  in  the 
counter-clockwise  direction,  the  change  in  0,  that  is  A0,  is  the  negative  of 
this  rotation  angle  6  (A0=  -5)  .  This  is  an  orthogonal  transformation,  and 
the  length  of  the  rotated  vector,  R',  is  the  same  as  the  length  of  the  original 
vector,  R. 

For  very  small  rotation  angles  sin  (5)  «  5 ,  and  cos  (5)  =  1  . 

Plugging  in  these  approximations  and  reversing  the  order  of  the  terms  in 
equation  (2),  we  have: 

(3)  X'  «  X  +  y5 

(4)  y’  =  y  -  x5 

Equations  (3)  and  (4),  along  with  a  third  equation  which  keeps  track  of 
the  cumulative  angle  of  rotation  (when  this  is  relevant),  are  the  main 
working  equations  of  the  CORDIC  algorithm.  The  details  of  this  procedure  are 
discussed  below  in  the  ALGORITHM  section  (Section  III). 

The  transformation  equations  are  now  no  longer  orthogonal,  and 
correspond  not  only  to  a  rotation,  but  also  a  stretching  of  the  vector  It  is 

shown  below  that  the  stretch  factor  (K)  equals  Vl  +5^ 


(5)  R'  *V{xf +(yf 

(6)  R'  =  V(x  +  y  5f  +  (y  -  X  5f 

(7)  R'  _  V x^  +  y25^  +  y2  +  x-5^ 

(8)  R-  =  V(x2+  y2)(l  +5^) 
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(9)  R’  =  R  V 1  +  5^ 

Furthermore,  5  no  longer  represents  the  angle  of  rotation  for  the 
vector,  but  instead  the  vector  will  have  been  rotated  clockwise  through  an 
angle  a  equal  to  the  arctan{5).  The  fact  that  a  equals  the  arctan(5)  is  proven 
next. 

Define  a  vector  V  that  has  the  same  length  as  R  and  the  same  direction 

as  R'. 


Figure  2:  Gordie  Rotation 


Since  the  magnitude  of  V  =  R  =  R'  /  Vl  +  6^  ,the  components  of  V, 
namely  and  ,  are  equal  to  x*  /  Vl  +  6^  and  y'  /  V 1  +  5^ .  respectively. 

Since  x’  =  x  +  y*6  ,  we  have  Xy  =  {  x  +  y*5 )  /  Vl  +  5^  ,  and  therefore, 

(10)  Xv  =  X  /  V77?  +  y  *  6  /  VlT? 

The  V  vector  is  the  R  vector  after  an  orthogonal  clockwise  rotation  through 
an  angle  a,  the  transformation  equation  for  Xy  has  the  form 

(11)  Xy  =  X  *  cos  (a)  +  y  *  sin  (a) 

Comparing  equations  (10)  and  (11)  for  Xy  we  see  that 

(12)  sin  (5)  =  5  /  Vl  +  5^  and  (13)  cos  (a)  =  1  /  V 1  +  5^ 

Recall  that 
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(14)  tan  (a)  =  sin  (a)  /  cos  (a) 

Plugging  the  expressions  in  (12)  and  (13)  into  (14)  we  get  tan  (a)  =  5  or 
a  =  arctan(5).  The  same  result  can  be  obtained  by  an  analysis  of  the  y 
component  of  V. 


III.  THE  ALGORITHM 

There  are  two  modes  for  the  CORDIC  algorithm.  One  is  called  vectoring; 
the  other,  rotation.  The  vectoring  mode  will  be  explained  in  detail  since  our 

problem  arose  in  this  mode  when  we  tried  to  compute  Por  an 

explanation  of  how  the  rotation  mode  can  be  used  to  compute  such  functions 
as  the  sine  and  cosine,  the  reader  should  consult  one  the  references'* -2.7.8, 
The  vectoring  mode  is  useful  when  the  x  and  y  components  of  a  vector 

are  given  and  the  magnitude  and/or  the  arctan(y/x)  are  desired.  In 

this  mode  the  successive  CORDIC  rotations  are  carried  out  in  such  a  way  as 
to  eventually  "force  y  to  zero".  Each  iteration  corresponds  to  a _ 

nonorthogonal  rotation,  and  stretches  the  vector  by  a  factor  of  V 1  +  5? . 

This  stretch  factor  is  independent  of  the  direction  of  the  rotation.  The 
cumulative  stretch  factors  are  listed  in  Table  2.  After  y  has  been  forced  to 
zero  (i.e.  the  vector  has  been  rotated  to  align  with  the  +x  axis  ),  the 

magnitude,  ^  j5  obtained  by  dividing  the  value  in  the  x  variable  by 

the  cumulative  stretch  factor. 

To  compute  Vx^T^  working  equations  are: 

(15)  Xj^i  =  X|  +  yi  5| 

(16)  yi^i  =  yi  -  Xj  5i 

where  for  the  ith  iteration  5j  =  ±  (1/2)'  and  i  =  0,  1,  2,  3  .  .  . 

The  ±  sign  is  selected  by  checking  whether  y\  is  positive  or  negative. 

In  order  to  force  y  to  zero,  if  yj  is  positive,  then  5j  is  positive,  and  XjSjis 
subtracted  from  y;  (  N.B.  Xj  is  always  positive  ).  Conversely,  if  y,  is 
negative,  then  6j  is  chosen  to  be  negative  also. 

Multiplying  Xj  or  y;  by  5;  is  achieved  by  right  shifting  the  value.  For 

example,  if  i  equals  3  then  83  equals  (1/2)3  .  jhe  value  of  y3  83  is  then 

computed  by  simply  shifting  the  binary  value  of  y3  three  places  to  the  right. 
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In  the  Vx2Ty2  computation  there  is  no  need  to  keep  track  of  the 
cumulative  rotation  angle.  However,  If  the  arctan(y/x)  of  the  original 
vector  is  desired,  then  one  simply  sums  up  the  angles  of  rotation  {a\) 
produced  by  each  iteration  (recall,  ai  »  arc  tan  (5i)  ). 
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IV.  INTEGER  ARITHMETIC  PROBLEM 

The  RTX  processor  is  equipped  with  specialized  square  root 
instructions.  This  routine  will  take  the  square  root  of  any  positive  integer 
up  to  31  bits  long  (corresponding  to  the  decimal  range  of  zero  to 
2,147,483,647).  This  may  seem  like  a  large  range,  but  in  the  special  case 

where  x  equals  y  in  ,  the  maximum  value  for  x  is  only  32,767.  This 

is  not  adequate  for  our  purposes.  We  tried  using  a  63  bit  square  root 
algorithm,  but  CORDIC  executed  faster.  Using  CORDIC  we  can  extend  the 
range  of  the  input  values,  x  and  y,  to  30  bits. 

Unfortunately,  when  we  tested  our  CORDIC  square  root  function,  we 
came  across  the  difficulty  illustrated  in  the  following  example^ _ 

Suppose  X  equals  333  and  y  equals  444  .  We  can  expect  to 

yijid  555  since  this  is  a  3-4-5  triangle.  Below  we  present  a  table  of  x,  ,  y,  , 
Xj  5i ,  and  yjSj  after  each  iteration  as  determined  by  the  algorithm  discussed 
above. 


Table  1:  Example 


i 

Xi 

yi 

yi5i 

0 

333 

444 

333 

444 

1 

777 

1  1  1 

388 

55 

2 

832 

-277 

208 

-70 

3 

902 

-69 

112 

-9 

4 

91  1 

43 

56 

2 

5 

913 

-13 

28 

- 1 

6 

914 

1  5 

14 

0 

7 

914 

1 

7 

0 

8 

914 

-6 

3 

-  1 

9 

915 

-3 

1 

- 1 

1  0 

916 

-2 

0 

- 1 

1  1 

917 

-2 

0 

-  1 

1  2 

918 

-2 

0 

- 1 

1  3 

919 

-2 

0 

-  1 

1  4 

920 

-2 

0 

-  1 

1  5 

921 

-2 

0 

-  1 

The  reader  will  note  that  after  iteration  #7  the  value  of  y  is  closest 
to  zero.  If  the  value  of  x  after  iteration  #7  (  namely,  914  )  is  divided  by  the 
stretch  facter  (  see  table  2  )  of  1.6466932543,  and  then  rounded  to  an 


6 


N  ADC-91 067-50 


integer,  the  result  turns  out  to  be  the  correct  integer,  555.  However,  ^rter 
iteration  #9,  y  is  stuck  at  -2,  but  x  (  and  therefore  the  result  )  continues  to 
grow. 


Table  2:  Stretch  Factors 


Iteration  Number 

Stretch  Factor  (K) 

0 

1.4142135624 

1 

1.5811388301 

2 

1.6298006013 

3 

1.6424840658 

4 

1.6456889158 

5 

1.6464922787 

6 

1.6466932543 

7 

1.6467435066 

8 

1.6467560702 

9 

1.6467592111 

1  0 

1.6467599964 

1  1 

1.6467601927 

1  2 

1.6467602418 

1  3 

1.6467602540 

1  4 

1.6467602571 

1  5 

1.6467602579 

1  6 

1.6467602581 

1  7 

1.6467602581 

1  8 

1.6467602581 

1  9 

1.6467602581 

20 

1.6467602581 

2  1 

1.6467602581 

22 

1.6467602581 

23 

1.6467602581 

24 

1.6467602581 

25 

1.6467602581 

26 

1.6467602581 

27 

1.6467602581 

28 

1.6467602581 

29 

1.6467602581 

30 

1.6467602581 

31 

1.6467602581 
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V.  SOLUTIONS 

We  considered  several  ways  to  patch  the  algorithm.  Since  the  y  value 
could  not  always  be  forced  exactly  to  zero,  we  needed  another  condition  that 
would  reliably  halt  the  iterative  process  without  introducing  too  much  error 

in  the  result  (Vx^  +  y2j  considered  checking  for  small  rates  of  change  in 

X,  y,  x5,  or  y5.  We  decided  instead  to  check  whether  the  absolute  value  i&f  y 
was  less  than  some  predetermined  cutoff  value  as  our  halt  condition.  The 
values  in  Table  1  suggested  to  us  that  if  the  absolute  value  of  y  became  less 
than  three,  it  was  time  to  stop.  This  condition  was  tested  by  looping 
through  millions  of  combinations  of  integers  that  maintain  the  3-4-5 
proportionality  and  were  in  our  range  of  interest.  We  also  decided  to  test 
other  limits  for  y.  The  limit  for  y  was  incremented  from  0  to  127.  Table  3 
is  a  representative  selection  of  the  distribution  of  errors  as  a  function  of 
the  |y|  cutoff.  The  error  frequency  counts  were  truncated  to  32760  to  avoid 
overflow.  When  the  ly|  cutoff  was  less  than  three,  a  second  peak  in  the 
error  distribution  appears  between  10  and  14.  These  occurrences  resulted 
from  cases  which  were  never  halted  at  maturity.  The  drift  from  the  correct 
result  continued  until  the  DO  loop  was  completed  (32  iterations). 

Using  the  combinations  of  integers  that  maintain  the  3-4-5 
proportionality,  the  error  stayed  below  six  for  a  broad  range  of  ly|  cutoff 
values.  Eventually,  at  a  sufficiently  high  cutoff  (approximately  100)  the 
size  of  the  error  began  to  rise  due  to  premature  halting  of  the  algorithm. 
These  cases  involved  small  initial  values  of  x  and  y.  In  particular,  when  the 
initial  value  of  y  was  less  than  the  cutoff,  the  algorithm  halted  immediately 
and  returned  the  initial  value  of  x  as  its  result. 
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Table  3:  Error  Frequency  Vs.  Size  of  Error  and  Cutoff 


Error 

lykO* 

1 

2 

3 

28 

60 

80 

101 

0 

26 

665 

1810 

3078 

32760 

32760 

32760 

32760 

1 

2566 

16177 

32760 

32760 

32760 

32760 

32760 

32760 

2 

23370 

32760 

32760 

32760 

32760 

28259 

19375 

19150 

3 

32760 

32760 

32760 

32760 

8528 

6803 

5446 

4359 

4 

21159 

26561 

19078 

9693 

1446 

734 

574 

436 

5 

8043 

6709 

3809 

2385 

61 

4 

2 

22 

6 

3159 

1188 

272 

138 

1 

0 

0 

1 

7 

1190 

432 

5 

1 

0 

0 

0 

0 

8 

812 

125 

0 

0 

0 

0 

0 

0 

9 

15729 

3839 

0 

0 

0 

0 

0 

0 

10 

32760 

21286 

8 

0 

0 

0 

0 

0 

1  1 

32760 

16511 

1  1 

0 

0 

0 

0 

0 

1  2 

7576 

3258 

3 

0 

0 

0 

0 

0 

1  3 

1465 

510 

5 

0 

0 

0 

0 

0 

1  4 

412 

239 

37 

0 

0 

0 

0 

0 

1  5 

68 

3  1 

1 

0 

0 

0 

0 

0 

1  6 

2 

0 

0 

0 

0 

0 

0 

0 

1  7 

1 

0 

0 

0 

0 

0 

0 

0 

1  8 

0 

0 

0 

0 

0 

0 

0 

0 

1  9 

0 

0 

0 

0 

0 

0 

0 

0 

*  This  is  equivalent  to  the  standard  CORDIC  algorithm  (no  lyl  cutoff). 

Other  combinations  of  integers  were  also  tested.  For  example, 
integers  that  maintain  the  5-12-13  proportionality,  as  well  as  integers 
generated  randomly,  were  studied.  The  general  features  of  the  distribution 
of  errors  as  a  function  of  the  ly|  cutoff  remained  the  same;  however,  the 
region  where  the  errors  were  less  than  six  moved  around. 

The  function  which  we  finally  implemented  involves  a  hybrid 

approach  to  evaluating  Vx^  +  y2_  Whenever  the  input  values  of  both 
X  and  y  are  smaller  than  32768,  the  RTX  processor's  31  bit  square 
root  function  is  employed.  Otherwise,  CORDIC  with  a  lyl  cutoff  of  1(X) 
is  used.  This  combined  the  best  of  both  worlds.  The  built  in  routine 
was  very  fast,  but  could  not  handle  large  numbers;  whereas,  CORDIC 
produced  a  much  smaller  per  cent  error  for  large  numbers  than  it 
did  for  small  numbers.  Setting  the  lyl  cutoff  at  100  has  the 
advantage  of  providing  a  relatively  quick  exit  condition. 
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Furthermore,  very  little  is  lost  with  this  choice  of  cutoff  since  we 
only  use  CORDIC  for  large  values  of  x  and  y.  Suppose,  for  example, 
the  initial  values  of  x  and  y  are  30,000  and  40,000  respectively. 

Since  one  of  these  numbers  is  larger  than  32768  we  would  utilize 

CORDIC.  The  expected  result  for  Vx^T^  jg  50,000.  When  the  vector 
has  been  rotated  such  that  y  =  100,  the  value  of  x  is  then  49,999.9 
(ignoring  the  stretch  factor  for  the  sake  of  argument).  The  truncated 
value  of  49,999  is  only  one  less  than  the  correct  value  of  50,000. 


VI.  CONCLUSION 

While  the  CORDIC  algorithm  provides  a  simple  method  of  evaluation  for 
a  wide  variety  of  functions,  we  found  that  caution  is  necessary  in  certain 

circumstances.  In  particular,  when  integer  arithmetic  is  used  and 
is  evaluated  by  CORDIC,  significant  errors  sometimes  arise.  This  is 
especially  bothersome  for  small  initial  values  of  both  x  and  y.  One  way  to 
handle  this  problem  is  to  place  a  cutoff  condition  on  the  absolute  value  of  y. 
Usually,  a  built  in  square  root  function  is  available;  however,  its  range  may 
be  too  limited.  We  recommend  using  the  built  in  function  because  of  its 
speed  and  accuracy  whenever  it  is  possible,  and  using  CORDIC  with  a 
suitable  cutoff  on  the  absolute  value  of  y  to  extend  the  range. 
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